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Abstract 



There is increasing evidence that many km-sized bodies in the Solar System are piles of rubble bound together 
by gravity. We present results from a project to map the parameter space of collisions between km-sized 
spherical rubble piles. The results will assist in parameterization of collision outcomes for Solar System 
formation models and give insight into fragmentation scaling laws. We use a direct numerical method to 
evolve the positions and velocities of the rubble pile particles under the constraints of gravity and physical 
collisions. We test the dependence of the collision outcomes on impact parameter and speed, impactor spin, 
mass ratio, and coefficient of restitution. Speeds are kept low (< 10 m s^^, appropriate for dynamically 
cool systems such as the primordial disk during early planet formation) so that the maximum strain on the 
component material does not exceed the crushing strength. We compare our results with analytic estimates 
and hydrocode simulations. We find that net accretion dominates the outcome in slow head-on collisions 
while net erosion dominates for fast off-axis collisions. The dependence on impact parameter is almost equally 
as important as the dependence on impact speed. Off-axis collisions can result in fast-spinning elongated 
remnants or contact binaries while fast collisions result in smaller fragments overall. Clumping of debris 
escaping from the remnant can occur, leading to the formation of smaller rubble piles. In the cases we 
tested, less than 2% of the system mass ends up orbiting the remnant. Initial spin can reduce or enhance 
collision outcomes, depending on the relative orientation of the spin and orbital angular momenta. We derive 
a relationship between impact speed and angle for critical dispersal of mass in the system. We find that our 
rubble piles are relatively easy to disperse, even at low impact speed, suggesting that greater dissipation is 
required if rubble piles are the true progenitors of protoplanets. 



1 INTRODUCTION 



There is growing interest in understanding the dynamics of collisions between small bodies in the Solar 
System. Typically such collisions are divided into two regimes: those dominated by material strength and 
those dominated by self-gravity (Holsapple 1994). The transition from the strength to the gravity regime may 
occur at body sizes as small as 250 m for silicates (Love and Ahrens 1996). In this paper we present numerical 
results from simiilations of collisions in the gravity regime. Our experiments are primarily concerned with 
low-speed collisions between equal-mass, km-sized "rubble piles" , gravitationally bound aggregates of loose 
material. We believe that these experiments will shod light on the collisional dynamics of the protoplanetary 
disk when typical encounter speeds are comparable to the surface escape speed (about 1 m s~^ for km-sized 
planetesimals of 2 g cm~^ bulk density). 

1.1 Definitions 

We begin with definitions of terms frequently encountered in the context of binary collision experiments. 
Typically in the literature one impactor (the larger one) is stationary and is considered to be the target 
while the other (the smaller one) is moving and is called the projectile. In our experiments, the impactors 
are comparable in size and are both in motion, so we generally do not distinguish between a target and a 
projectile. Most laboratory experiments involve solid targets that possess tensile strength, so the outcome is 
measured in terms of the extent of disruption or shattering of the target. A critical or catastrophic shattering 
event is one in which the largest post-impact fragment (the remnant) has 50% of the target mass. Following 
a recently adopted convention in the literature (Durda et al. 1998), we use Qg to denote the kinetic energy 
per unit target mass to achieve critical shattering. A rubble pile, by definition, has no tensile strength, so 
Qg is effectively zero. However a rubble pile can still be disrupted or shattered in the sense that one or more 
of the component particles becomes separated from the rest for at least an instant. 

For collisions in free space, fragments or particles are said to be dispersed if they attain positive orbital 
energy with respect to the remnant. Hence a critical or catastrophic dispersal is one in which the largest 
remnant is left with 50% of the original target mass after the remaining material has dispersed to infinity. The 
energy per unit target mass to achieve this is denoted by Q^. In our experiments, since we do not distinguish 
between a target and a projectile, Q^j refers to the energy per unit total mass, in the center-of-mass frame, 
needed to critically disperse the entire system. 

Finally, we define erosion to mean permanent removal of mass from a body, and accretion to mean 
permanent retention of mass. In the context of our experiments, net erosion means that one body (the 
largest if the impactors are of unequal mass) had less mass at the end of the run than it started with. Net 
accretion means it had more mass at the end. 

1.2 Motivation 

Many asteroid characteristics are inconsistent with monolithic configurations. Recent observations by the 
Near Earth Asteroid Rendezvous spacecraft of Mathilde, a 53-km C-class asteroid, are particularly suggestive. 
First, Mathilde's largest crater is enormous: it has a diameter of 33.4 km, almost 7 km larger than the 



asteroid's mean radius (Veverka et al. 1997). Numerical hydrocode simulations and laboratory experiments 
strongly suggest that in order for Mathilde to have survived the impact that formed such a substantial 
crater, the asteroid must be made of some material that does not efficiently transmit energy throughout the 
body (Housen and Holsapple 1999). Second, Mathilde has a remarkably low density of 1.3 g cm^-^ (Yeomans 
et al. 1997), about one third the average value for the chondritic meteorites that are thought to originate 
from C-class asteroids (Wasson 1985). Such a low density suggests that Mathilde is highly porous. If true, 
the voids in the material could impede the transmission of energy from a collisional shock wave and allow 
a rather weak body to survive an otherwise catastrophic impact event (Asphaug et al. 1998, Housen and 
Holsapple 1999). 

In addition to Mathilde, the surfaces of 243 Ida, 951 Gaspra, and Phobos show several sizable craters 
that have diameters on the order of the mean radius of the body (for references, see Richardson et al. 1998, 
hereafter Paper I). As in the case of Mathilde, the energy necessary to create craters of this size would 
disperse or disrupt the original body if it were solid (Asphaug et al. 1998). 

Further evidence for the prevalence of rubble piles comes from asteroid spins. In a sample of 107 asteroids 
smaller than 10 km in diameter, Harris (1996) found that the spin period distribution truncates at fast spin 
rates, where rubble piles would start to fly apartQ 

One explanation for the observed characteristics of these asteroids and their craters is that they are rubble 
piles. Although rubble-pile configurations are more susceptible to disruption by tidal forces than monolithic 
configurations (Paper I) , there is increasing evidence that rubble piles have a higher impact strength (Ryan 
et al. 1991, Love and Ahrens 1996, Asphaug et al. 1998). There are two scenarios for creating a rubble-pile 
asteroid: 1) the asteroid is initially one solid body of material and is rubblized over time by multiple impacts; 
2) the rubble-pile configuration of the asteroid is primordial. Regardless of how rubble-pile asteroids are 
formed it is interesting to investigate how they interact and evolve in the Solar System. 

1.3 Laboratory Experiments: Strength vs. Gravity 

Ryan et al. (1991) presented results from a laboratory study of impacts into weak inhomogeneous targets. 
Due to practical limitations they used ~ 0.5-cm targets of gravel and glue. As a result, their specific 
experimental results are firmly rooted in the strength regime. However, the most general conclusion that 
the group arrived at from dropping, crashing, and shooting at the gravel aggregates was that the relatively 
weak targets have a surprisingly high impact strength. In other words, it took a large amount of energy (at 
least Q*g — 2.6 x 10^ J m"'^) to critically disrupt or shatter the target such that the largest remnant was one 
half the mass of the original object. The nonuniformity of the target causes a greater fraction of the impact 
energy to dissipate thermally; therefore, the collisional shock wave is more efficiently absorbed by the target. 

Laboratory experiments on Earth to investigate directly the collisional dynamics of the gravity regime 
are difficult to conduct since the target size necessary to reach this regime is impractically large. Instead, 
overpressure and centrifuge techniques have been used to artificially simulate the gravity regime in the 
laboratory. In an overpressure experiment, Housen et al. (1991) used nitrogen gas at various pressures to 

^At least one asteroid spinning faster than this hmit has since been discovered (Ostro et al. 1999), but its small size (~ 30 
m) puts it comfortably in the strength regime. 



mimic the lithostatic stress felt inside a large target. At these pressures they were unable to carry out true 
impact tests so they used a buried charge instead of a projectile. As the pressure was increased, the size of 
the largest remnant after each explosion also increased, indicating a transition from the strength-dominated 
regime to the pressure-dominated regime. Housen et al. (1991) argued that the pressure regime was analogous 
to the gravity regime and extrapolated a scaling law for the gravity regime from the overpressure data. This 
laboratory study has two important drawbacks: 1) by using a buried charge the experiment does not model 
the actual surface dynamics of an asteroid during an impact; 2) the gas overpressure is not an r~^ force law. 
They were able to reach a regime in the laboratory that was not dominated by the strength of the material, 
but it is unclear whether the gravity-regime scaling law derived from the overpressure data is valid. 

In a centrifuge experiment, Housen and Holsapple (1999) were able to conduct true impact tests by firing 
a small projectile (a polyethylene cylinder 0.65 cm in diameter) from a gas gun strapped to the arm of a 
centrifuge. They positioned a porous target (composed of quartz sand, perlite, fly ash, and water) at the 
end of the arm. The centrifuge was used to mimic the gravitational force at the surface of a much larger 
body. The use of the centrifuge introduces second-order complexities due to the Coriolis force and the field 
orientation in general at the surface of the cylindrical target (though this is only really a problem in the event 
of high ejecta trajectories). In addition, the fiat surface of the target may subtly affect crater morphology. 
Nonetheless, this experiment showed that porous targets in the gravity regime are efficient at absorbing 
impact energy at the surface by compacting the underlying material. 

1.4 Numerical Simulations of Collisions and Tidal Disruptions 

Extrapolations of laboratory experiments have resulted in rough strength and gravity scaling laws. In 
order to truly understand the coUisional dynamics and evolution of large bodies, numerical simulations are a 
necessity. Love and Ahrens (1996) used a three-dimensional "smoothed particle hydrodynamics" (SPH) code 
to simulate high-speed catastrophic collisions. They used various impact speeds (3 7 km s^^), impact angles 
(5-75°), target diameters (10-1000 km), and projectile diameters (0.8-460 km) in order to explore a large 
region of parameter space. The big targets placed the experiments securely in the gravity regime, allowing 
the researchers to treat gravity carefully and neglect the strength and fracturing of the target completely. 
Their extrapolated scaling law for the gravity regime placed the transition from the strength to the gravity 
regime at a target diameter of 250 ± 150 meters, much smaller than that predicted by laboratory experiments 
(Holsapple 1994). Love and Ahrens (1996) argue that since smaller asteroids are more common than larger 
ones, a given asteroid is more likely to suffer a shattering impact before a (lisp(!rsing impact. Thus, it seems 
plausible that many asteroids in our Solar System are at least partial rubble piles. 

More recent simulations have had similar results. Asphaug et al. (1998) conducted three high-speed (5 
km s~^) impact experiments using a solid target, a partially rubblized contact binary, and a totally rubblized 
ellipsoidal target. In each case the researchers used a small projectile six orders of magnitude less massive 
than the target. There are three major conclusions from this study: 1) it is much easier to disrupt a solid 
target than it is to disperse it — this conclusion is evidence that it is possible to change a solid body into a 
rubble pile with impacts; 2) rubble regions can insulate and block energy from traveling through a body — in 
a contact binary, for example, one end could be critically disrupted while the other remains undamaged; 



3) the fully rubblized targets efficiently localize the energy transmitted during a collision which in turn 
minimizes the damage outside the collision region and allows weak bodies to survive high-energy impacts 
with much less damage than solid targets. This again implies that many small bodies in the Solar System 
may be rubble piles. 

Watanabe and Miyama (1992) used 3-D SPH code to investigate the effects of tidal distortion and shock 
compression from collisional impacts in the process of planetary accumulation. They used two equal-sized 
spherical bodies and assumed a perfect Newtonian fluid. It is important to note that their code did not 
model an incompressible fluid (their adopted polytropic indices were always greater than zero). As a result 
of experimenting with impact angle, speed, and density gradients, they found that tidal forces can enlarge 
the coalescence rate of planetesimals by almost a factor of 2. In addition, when the initial speed of the 
impactor is signiflcantly lower than the escape speed of the system, less than a few percent of the total mass 
is lost from the system in the collision. They did not attempt any simulations with initial speeds in excess 
of 50% of the escape speed. 

In Paper I, Richardson et al. numerically simulated the effects of the Earth's tidal force on rubble-pile 
asteroids. Unlike Watanabe and Miyama (1992), they simulated the Earth-crossing asteroids as incompress- 
ible fluids using a hard-sphere model. They varied the asteroids' speed, spin, shape, and close- approach 
distance. Generally, slow-moving, close-approaching, prograde-rotating, elongated asteroids were the most 
susceptible to tidal disruption. They found several distinct classes of outcome: in the most violent disruption 
cases, the asteroid was stretched into a line and recollapsed into a "string of pearls" reminiscent of Comet 
D/Shoemaker-Levy 9 at Jupiter; for moderate disruptions, large pieces of the asteroid were stripped off in 
many cases, forming satellites or contact binaries; the mildest disruptions resulted in little mass loss but 
significant shape changes. These various outcomes could lead to the formation of crater chains (Bottke et al. 
1997), asteroid satellites and doublet craters (Bottke and Melosh 1996a,b), and unusually shaped asteroids 
(Bottke et al. 1999). 

Durda (1996) carried out simulations to study how readily satellites form as a result of mutual gravi- 
tational attraction after the catastrophic disruption of the progenitor. Durda (1996) came to three major 
conclusions: 1) satellites do form immediately after a catastrophic collision; 2) contact binaries form more 
easily than true binary systems; 3) the binary systems form in a wide range of size ratios. It is important 
to realize that Durda (1996) assumed a power-law mass distribution for the catastrophically fragmented 
asteroid. The slope index used (1.833) was taken from extrapolations of laboratory experiments. 

1.5 Implications for Planet Formation 

Traditionally, numerical simulations of planet formation use extrapolations of impact experiments in the 
strength regime to model the effects of fragmentation in planetesimal collisions (e.g., Greenbcrg et al. 1978; 
Beauge and Aarseth 1990; Wetherill and Stewart 1993). From what we have already seen, such extrapolations 
may give misleading results since generally much more energy is needed to disperse than to disrupt a 
planetesimal in the gravity regime. Moreover, effects of impact angle, spin, and impactor mass ratio are 
often not taken into account. In the case of rubble piles, no empirical model actually exists. For example, we 
might expect reaccumulation like that seen in the tidal disruption models to also occur after the catastrophic 



impact of two rubble-pile planetesimals. In this paper we aim to explore these issues by simulating collisions 
between rubble-pile bodies over a wide range of parameter space and determining the implications of the 
results for planet formation. In Section |21 we describe our numerical method and analysis technique. Our 
results are presented in Sectional followed by a general discussion in Section 0] We give our conclusions in 
Section 

2 METHOD 

The simulation and analysis of the collisions presented here combine numerical methods introduced in Paper 
I and Richardson et al. (1999), hereafter Paper II. The rubble pile model is an extension of the model used 
for studying the tidal disruption of asteroids (Paper I) . The integration engine is an extension of the parallel 
tree code used for planetesimal evolution simulations (Paper II). 

2.1 Rubble Pile Model 

Each rubble pile in our simulations consists, at least initially, of a fixed number of equal-size hard spheres 
arranged in "hexagonal close-packed" (HCP) form. The rubble piles are typically generated by specifying the 
bulk semi-axes, bulk density, and approximate number of particles (alternatively, the particle radius and/or 
density can be used as independent parameters) . The generator attempts to match the requested properties 
on the basis of the estimated HCP efficiency of a sphere as a function of bulk radius or number of particles 
(derived from power-law fits to our own numerical experiments). Once the rubble pile is constructed, the 
constituent particles are reduced in size by a fixed factor (usually 1%) and given a small random velocity 
kick (no more than 10% of the particle surface escape speed in magnitude). This is to facilitate attaining the 
initial equilibrium (c/. Section l2.4|l . Finally, the rubble pile is tagged with a unique "color" so that mixing 
can be studied visually and statistically. 

The collisional properties of the constituent particles are specified prior to each simulation. These include 
the normal and tangential coefhcients of restitution, e„ and et (cf. Richardson 1994). Except for certain 
explicit test models, these values generally were fixed at e„ = 0.8 (mostly elastic collisions with some 
dissipation) and et — 1.0 (no surface friction). Bouncing was the only possible collision outcome: no mergers 
or fragmentations of particles were allowed. The value of e„ was chosen to be consistent with Paper I and 
is similar to experimentally determined values used in the literature (e.g., Beauge and Aarseth 1990). Note 
that in the perfectly elastic case, particles cannot re-collapse into condensed rubble piles after a disruption 
event but instead completely disperse or at best end up in centrally concentrated swarms. In the case of 
tidal disruption (Paper I) the outcome is relatively insensitive to the choice of e„, so long as e„ < 1. For 
the present study, however, varying e„ has a stronger effect, an issue we explore in Section 13.21 We did not 
include surface friction in the present study, in order to keep the number of test cases manageable. 

There are two circumstances under which e„ is allowed to change. First, if the relative speed of two 
colliding particles is less than 10% of their mutual escape speed (i.e., typically ~ 1 cm s^^), e„ is set to unity 
(no dissipation). This is to prevent computationally expensive "sliding motions" (Petit and Hcnon 1987). 
Second, if the coUision speed exceeds 10 m s^^, e„ is set to 0.2 (highly dissipative) . This is to crudely model 



damping through internal fracture as the impact stress pvc {p = internal density, c = sound speed ^ IQ^ 
m s^^) exceeds the "rock" strength 10^ N m^^). This is not intended to be a physically rigorous model 
but rather a simple mechanism to prevent unrealistically high collision speeds. Initial encounter speeds 
between rubble piles were generally kept closer to 1 m s^^ in any case. Also, particle sizes were kept roughly 
comparable across rubble piles in order to minimize any strength- uersws-size biases. 

It is important to note that neither rolling nor true sliding motions are modeled in our code. Moreover, 
particles cannot remain mutually at rest in contact (i.e., there are no surface normal forces). Instead, the 
constituent particles of an otherwise quiescent rubble pile are in a constant state of low-energy collisional 
vibration (dictated by the minimum sliding condition described above). Nevertheless, such small bounces 
can mimic transverse motions in an approximate sense in the presence of shear flow, giving realistic bulk 
properties to the material. To test this we have simulated the formation of sand piles using our collision code 
(with surface friction) that give reasonable values for the angle of repose when compared with laboratory 
experiments. 

2.2 Numerical Code 

Our simulations were performed using a modified version of a cosmological TV-body code, pkdgrav (Stadel 
and Quinn, in preparation; data structures described in Anderson 1993, 1996). This is a scalable, parallel tree 
code designed for ease of portability and extensibility. For the parameter space study, the parallel capability 
was not exploited owing to the modest number of particles in each run (a few thousand). However, even 
in serial mode, pkdgrav is arguably more efficient than any other existing code with similar capability. In 
particular, it is superior to box_tree, the code used in Paper I, which could only handle a few hundred 
particles in practical fashion. 

A low-order leapfrog scheme is used as the pkdgrav integrator. The comparative simplicity of this scheme 
is a big advantage for collision prediction since particle position updates are linear in the velocity term. This 
means that every possible collision within the timestep can be determined in advance and in the correct 
sequence. Timesteps are smaller than in higher-order schemes for the same accuracy, but the cost of each 
gravity calculation is far outweighed by the collision search once the rubble piles are in contact, and is 
comparable otherwise. Moreover, away from collision, particle trajectories are integrated symplectically, 
eliminating spurious numerical dissipation. For further detail and references, refer to Paper II. 

Although the collision search is relatively expensive, the scaling is modest: O(A^logA^) with particle 
number and linear with the number of collisions per interval. A typical encounter between thousand-particle 
rubble piles can generate ~ 10^ collisions over the course of a run! A balanced k-d tree (Bentley and Fried- 
man 1979) is used to search for possible collisions at the beginning of each timestep, giving the 0{N \ogN) 
dependence. Once a collision is performed, only particles that might be affected by the event in the remaining 
interval (numbering typically ^ N) are reconsidered via the neighbor search, giving the near linear depen- 
dence on the number of collisions. This latter enhancement is an improvement to the Paper II code, which 
did not require as much sophistication given the low collision frequency per step. Note that the collision 
search can also be performed in parallel, which proved necessary for the large- TV models presented in Section 
13.31 below. 



2.3 Hardware 

The parameter space models were lun on a local cluster of 16 300-MHz Intel Pentium lis using the High 
Throughput Computing (HTC) environment condor (c/. http://www.cs.wisc.edu/condor/), version 6.0, 
under RedHat Linux 5.0 with a 2.0.35 Linux kernel. The condor system supports automatic scheduling, 
submission, and restarting of jobs on shared resources, greatly simplifying management. A typical run 
required between 12 and 72 wallclock hours to complete and each generated ~ 25-50 MB of data. Models 
requiring parallel resources were run either on a local cluster of 4 433-MHz DEC Alpha PCs connected with 
a fast ethernet switch, or on a local SGI Origin 200 with 4 180-MHz processors running IRIX 6.4. Both 
platforms typically achieved sustained performances of several hundred megaflops. 

2.4 Initial Conditions 

Generation of initial conditions and analysis of results were performed using code auxiliary to pkdgrav. The 
rubble pile generator has already been described fSection l2.1|l . Each new rubble pile was first run in isolation 
(with or without spin) using pkdgrav until the velocity dispersion of the constituent particles achieved a 
stable equilibrium. Next a new "world" was created by using a small program to position and orient any 
number of equilibrated rubble piles (always two in the present study) prior to simulation. Spherical bodies 
were usually given a random orientation in order to reduce the effect of HCP planes of symmetry. Bulk 
velocities were then applied to each rubble pile. Other rubble pile properties that could be changed at this 
point included the total mass, bulk radius, bulk density, and color. For the exploration of parameter space, 
usually only the positions (in the form of y offsets), velocities, spins, and colors were modified. Once all the 
rubble piles were in place, the world would be adjusted so that the center of mass coincided with the origin 
and the velocity of the center of mass was zero. The output world was then read in by pkdgrav and the 
simulation would begin. 

To facilitate the exploration of parameter space, a series of Unix scripts were written to generate and 
monitor each run. Starting with a given pair of rubble piles and a list of desired initial impact parameters, 
speeds, and spins, the world generator was run automatically to create the necessary initial conditions and 
support files in separate run directories. The scheduler condor was then invoked to farm the jobs out to all 
available machines. Analysis was performed on the fly using a machine outside the condor pool for maximum 
efficiency. 

The choice of initial conditions was governed largely by prior test simulations. For the parameter space 
exploration, 10 values of impact parameter b and 10 values of initial relative speed v were chosen for each set 
of runs, where a set consisted of a fixed choice of spin and/or offset direction (see Section 13 for a complete 
description of each model). From the test simulations it was clear that only about half of the possible 100 
runs for each model were needed to find the representative cases and the Q*jj boundary. In a plot of b vs. v, 
the important region is the lower- left triangle (see Fig. El for an example). The b and v values were therefore 
chosen to sample this region as finely as possible in a practical amount of time. Models with spin were chosen 
to sample representative combinations of spin and orbital angular momentum at a fixed rotation period. 



2.5 Coordinate System and Units 

We use an inertial Cartesian coordinate system in free space for our simulations, with the origin at the center 
of mass. In the parameter space studies, the initial motion of the colliding bodies is in the ±x direction. 
Any initial impact parameter is measured in the ±y direction. Most debris actually travels in directions 
perpendicular to the original axis of motion (c/. Section ^31- 

A natural unit for the impact parameter b is the sum of the radii Ri + R2 of the two (spherical) impactors. 
Hence 6 = implies a head-on collision while 6 = 1 is a grazing encounter. Note however the true trajectories 
will generally be hyperbolae; no allowance is made for this in the definition. Since tidal effects may play 
an unpredictable role anyway, we adopt the simpler definition. In the absence of trajectory deflection, the 
impact angle is then </> — sin~^ 6, for b < 1. 

The unit for the initial relative speed v is more complicated. We chose a system in which v — indicates 
no relative motion and w = 1 is the estimated critical speed for dispersal. The critical speed is found by 
equating the initial total kinetic energy with the gravitational binding energy of a rubble pile made up of a 
spherical and homogeneous mixture of both colliders: 
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v...-M^ — , (1) 

where M is the combined mass, G is the gravitational constant, ^ is the reduced mass M1M2/M, and R is 
the radius of the sphere that contains the combined mass, assuming the same bulk density: 

R={Rl + Riy^\ (2) 

Note that the actual speed at impact will slightly exceed v due to gravitational acceleration. 

In the parameter space models, the initial separation in x for all cases was ~ 6R, effectively 2.5 Roche 
radii for the combined mass, i.e., far enough apart that initial tidal effects were negligible. The total energy 
of the system was positive in all cases. For completeness, the speed at infinity is given by 
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and the speed at impact is 



2.6 Run Parameters 
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Most pkdgrav run parameters assumed default values for these simulations (c/. Paper II). However, in 
addition to the collision parameters described in Section l2.ll the run time, timestep, and output frequency 
were specified explicitly for each model. 

The run time was initially 10 times the characteristic time 



where x is the initial separation. In most cases this is sufficient time for the post-collision system to reach a 
steady state. Some cases were run longer (typically a factor of 2) if necessary, on the basis of visual inspection 
of animations. 



The timestep for each run was set to a smaU vahie tg times a heuristic scale factor of 1/ {2v +1), arrived 
at by trial and error from our test runs (recall v is the initial speed, so to is a simple constant). The scaling 
ensures finer intervals for neighbor searches in higher-speed impacts (this is necessary to avoid missing any 
potential collisions). For our runs, to = 10^^ yr/27r, or roughly 50 s. Note that for objects with bulk density 
a few g cm^"^ the dynamical time 1/^/Gp ~ 1 h, comfortably large compared to the maximum adopted 
timestep. Generally our simulations are limited by the time needed to deal with particle collisions, so the 
gravity calculations can be of higher accuracy with little additional cost. 

Finally, the output frequency was chosen so that there would be about 200 outputs per run, suitable for 
smooth animations and analysis. 

2.7 Analysis Method 

Much of our analysis method is similar to that presented in Paper I; the reader is referred to that work for 
additional details. The basic strategy is to identify the largest post-collision remnant, compute its various 
properties, and generate statistics for the relative distribution of the smaller fragments. We use a slightly 
different clump- finding algorithm fSection l2.7.1(l . and now employ a shape drawing technique fSection l2.7.2|l . 
We have made other refinements that should improve the accuracy of the analysis. 

2.7.1 Clump finding 

The clump-finding algorithm iteratively refines guesses as to what constitutes a rubble pile by merging groups 
of particles together in bottom-up fashion. The first guess is that every particle in the system is its own 
rubble pile. On each pass basic properties are computed for each clump: mass, position, axis lengths, and 
orientation. Clumps are then compared in pairwise fashion. In order for two clumps to be merged (i.e., to be 
considered one clump), spheres of diameter equal to the major axes times a fixed linking scale and centered 
on each clump must overlap. If the scaled minor axes also overlap, then the clumps are merged. Failing 
that, if either body has its center of mass in the other's scaled ellipsoid, the bodies are merged. Otherwise 
no merge occurs. This process iterates until there are no more mergers during an iteration. 

This method is purely geometrical: gravitational groupings are not considered. This was done mostly 
because there is no natural gravitational length scale in the present context, unlike in Paper I where the Hill 
radius could be used. However, osculating elements of groups measured with respect to the largest fragment 
are still calculated and give a good indication of the future evolution of the system. The present method also 
differs from Paper I in that it treats each clump as an ellipsoid rather than a sphere, allowing more refined 
boundaries to be drawn. Through trial and error (visual inspection) a linking scale of 1.1 was adopted. 

2.7.2 Shape drawing 

During the course of the present investigation we came across some unusual, often asymmetrical shapes 
following collision events. In order to characterize these forms, a shape-drawing algorithm was devised. The 
algorithm attempts to trace the outer surface of a given rubble pile (either in cross section or by projection 
to the x-y plane). The resulting shape is equivalent to what would be measured by laser beams aimed at the 
surface in the direction of the center of mass. Note that this means any outcroppings can conceal underlying 



structure. Generally such complex surfaces are not seen in our models however (as confirmed by 3-D VRML 
viewing). The projection method is used in the parameter space plots of Sectional 



2.7.3 Mixing 

The unique color assigned to each rubble pile makes it easy to assess visually the degree of mixing following 
a collision. In order to make a more quantitative assessment, we have constructed the following statistic: 

1/2 



/mix 1 



E 



^c, world 



(6) 



^ X-^c' -v ^^c' ^c', world 

where subscript c denotes a color, subscript v denotes a subvolume of the rubble pile, Ny is the number 
of subvolumes, and "world" refers to the entire population of particles in the system. Note that particle 
number is conserved so that Ec 'tIc, world = M, the total mass of the system. This formula is generalized 
for any number of components (colors); in the present study only two populations are considered. A mixing 
fraction of unity implies the rubble pile contains a perfectly homogeneous mixture of the world colors. A 
value of zero means no mixing has taken place at all. 

Spherical subvolumes are used to sample different regions of the rubble pile (which itself need not be 
spherical). The size of the sample region is set so that it contains \/N particles on average. The center of 
a subvolume is chosen randomly within a rectangular prism enclosing the rubble pile. A new subvolume is 
chosen if the region is found to contain fewer than N^^^ particles. Otherwise, the argument of the Ei, 
Eq. is computed and added to the running sum. This is repeated until ^/N subvolumes are successfully 
sampled. 



3 RESULTS 

We now present the results of our simulations. First we describe the parameter space exploration which 
consisted of numerous runs of modest size. Highlights are shown in Fig. ^ where we have endeavored to 
illustrate the various classes of outcomes. Second we show the dependence on the coefficient of restitution 
e„ for a particular high-energy run. Finally we present the results of two high-resolution cases and compare 
with the corresponding moderate-resolution runs. 

3.1 Parameter Space 

We divided our exploration of parameter space into three models: a generic case as a baseline, a case with 
spinning impactors, and a case with unequal-mass impactors. Graphical summaries of these models are given 
in Figs. [3 and 01 which are discussed in detail below. 

3.1.1 Model A: Equal size, no spin 

Model A, our generic case, consisted of two equal-size rubble piles of 1 km radius and 2 g cm~'^ bulk density. 
The rubble piles were generated and equilibrated using the process described in Section 12.41 Each rubble 
pile contained 955 identical spherical particles of 83 m radius, so the packing efficiency was 55%^ The 

■^The effective packing efficiency is less than the maximum close-packed efficiency of 74% due to finite-size effects (Paper t). 



parameter space runs from 0.00-1.25 in b and 0.52-2.50 in v (the units of b and v are defined in Section 
12.51 Wcrit ~ 2.06 m s^^ for this model). The impact parameter values were chosen to encompass a range of 
dynamic interactions from head-on collisions to glancing distortions. The lowest value of v is twice the value 
corresponding to Wqo — (c/. Eq. smaller v leads to strong trajectory deflections). The largest value 
of V was chosen to be a moderately high-speed impact to ensure that the catastrophic dispersal regime was 
entered. 

Figure [5] summarizes the results of this model (Figs, nia)-(c) give snapshots of three distinct outcomes). 
The shapes in Fig. |21 trace the projected silhouettes of the largest post-encounter fragment at the end of 
each run. We have used nested squares of different line styles to divide our results into three mass regimes. 
A solid-line inner square indicates that the largest fragment contains 90% or more of the total mass of the 
system, i.e., nearly perfect accretion. A dashed- line inner square indicates that the largest fragment contains 
at least 50% but less than 90% of the total mass. The remaining cases contain less than 50% of the total 
mass in the largest fragment, i.e., net erosion. Note if there is no mass loss or exchange during the encounter 
the largest fragment will contain 50% of the total mass of the system by definition. We see in this model 
that 18 out of 55 runs (33%) result in net mass loss, although we caution that several cases are just on the 
border of 50%. 

The general trends in Fig. [21 are twofold, namely, as the encounter speed increases, the size of the largest 
fragment decreases, and as the impact parameter increases, the axis ratio increases, up to a certain point. 
Higher encounter speeds imply larger kinetic energy so it is more likely for the system to become unbound. 
Larger impact parameters imply larger angular momentum which results in an increase in the axis ratio 
until the critical spin value of the combined rubble pile is reached. In addition to the general trends, the 
middle-mass group has two distinct populations that reflect their formation history. The small 6, large v 
group (top left in the figure) represents a net loss of mass of 10%-50% from the system. The large b, small 
V group (lower right) represents grazing collisions with little mass loss or exchange. 

We note that for the head-on case our definition of Wcrit does not correspond to critical dispersal, rather, 
critical dispersal seems to occur at ^ 1.9wcrit- This probably reflects the fact that the energy of the collision 
is not immediately transported to all of the particles and that the voids in between the particles decrease 
the efficiency of energy propagation. Moreover, we did not take into account e,i in the definition of Wcrit- 
Regardless, Vcnt is intended as an approximate scaling only. 

More detailed results for this model are given in Table HI In the table, b and v have the usual definitions; 
Mi-om is the mass fraction of the largest post-encounter remnant; P is its instantaneous spin period in hours; 
e is the remnant's "ellipticity" : e = 1 — ^{q2 + Qs), where q2 = 02/01, = a^/ai, and ai > 02 > 03 are the 
semi- axes (e = is a sphere); /mix is given by Eq. JHJ; Mace, Morb, and Afcsc are the mass fractions that are 
accreting, orbiting, and escaping from the largest remnant, respectively!^ and ni, 77.2, and n are the number 
of single particles, two-particle groups, and discrete rubble piles (i.e., groups with three or more particles), 
respectively, at the end of the run. The Afrom column of Table compliments Fig. [5] by providing a finer 
gradation of the remnant mass. Note that Micm + M^cc + A/orb + Mcsc — 1- 

^To be considered accreting, a clump must have q < r + R, where q is the close-approach distance to the remnant, and r 
and R are the radii of minimal spheres enclosing the clump and remnant, respectively. This differs somewhat from Paper I. 



Table ^ shows how the remnant spin P is coupled to the angular momentum in each run. Since there 
are no external torques in the system, angular momentum is conserved. In the case of head-on collisions 
(6 = 0), there is exactly zero total angular momentum, which accounts for the large remnant P values (i.e., 
low spin). P is never infinite in these cases because some particles escape and carry angular momentum 
away from the remnant, even in the slowest collision case (w — 0.52). At higher collision speeds, more mass 
is carried away from the system, generally resulting in smaller P values. As b increases, so does the net 
angular momentum, resulting in faster spins (smaller P). This trend continues until 6^1 which corresponds 
to a grazing collision. In this case, the encounter generally does not result in a merger so the "remnant" is 
effectively one of the initial bodies plus or minus some mass exchange. Mass exchange and/or tidal torquing 
following deformation converts orbital angular momentum into spin angular momentum. As b increases 
further, there is little spin-up, since torquing becomes less effective. All of these trends can be seen in the 
table. 

Similarly, e depends on the total angular momentum of the system. Larger angular momentum allows 
the remnant to support a more elongated shape as long as most of the system mass ends up in the remnant 

75%, from the table). Consequently there is also a relationship between e and P: smaller P values 
correspond to larger e values, in general. The smallest P in the table is 4.1 h with e = 0.26; the largest e is 
0.45 with P — 4.3 h. These values are within the classical limit for mass retention at the surface: 



where p is the bulk density and we have assumed a2 ~ a^j. In this expression Pcrit — 2.3 h for a spherical 
rubble pile with p = 2 g cm"'^, and increases to infinity as e ^ 1. 

The sixth column in Tabled labeled /mix, gives the mean percent mixing fraction and standard deviation 
after 100 repeated measurements (recall the mixing calculation subvolumes are chosen randomly — cf. Section 
I2.7.3|l . The errors are a small fraction of the mean except when the mixing fraction itself is small. In the 
head-on case, /mix shows a simple trend of generally increasing with impact energy with a dip at the highest 
energy probably due to increased statistical fluctuation (the remnants are smaller). For most b values the 
situation is more complicated depending on whether the impactors accrete into a single body or exchange 
mass while remaining two separate bodies. For the largest &, little mass is exchanged, so the bodies remain 
essentially unmixed. 

The next three columns give dynamical information about the remaining mass of the system, i.e., the 
material not incorporated in the largest remnant. Generally most of this mass is escaping from the largest 
remnant (A/csc)- Typically only small amounts (< 10%) of mass are accreting (Mace) and/or orbiting (A/orb)- 
In two cases, however, Afacc is close to 50%; these are instances of near escape that were too computationally 
expensive to run until final accretion and represent the transition from a high- to medium-mass remnant. 

The final three columns contain information about the particle groupings at the end of each run. The 
number of free particles (ni) increases dramatically with v, but decreases with b. This trend is also seen in 
the number of two-particle groups (n2 ) and discrete rubble piles (n) . Groups can form either from accretion 
among the free particles due to gravitational instability or from being stripped off as a clump during the 
collision event. Note that n is always at least 1 because the remnant is included. 




(7) 



In summary, the outcomes of this model depend in a natural way on the total angular momentum and 
the impact energy of the system (both related to b and v). Larger b results in more elongated remnants 
with higher spins and reduced mixing. Larger v results in greater mass loss and increased mixing. In the 
remaining sections we explore how these trends are modified for non-identical or spinning bodies. 

3.1.2 Model B: Equal size, spin 

In Model B we added a spin component to the impactors. The spin vectors are oriented perpendicular to 
the orbital plane (i.e., along the ±z axis). The rotation period of the impactors is 6 h, the median rotation 
period of Near Earth Asteroids (Bottke et al. 1997). We investigated three cases: in Model Bl the spins of 
the impactors have opposite orientation; in Model B2 and B3 the spins have the same orientation but the 
impactors have opposite y offsets (Fig. Ol . By symmetry, these cases test all the unique z angular momentum 
combinations (spin + orbital). The remaining parameters are identical to those in Model A. 

Figures 01a)-(c) summarize the results for Models Bl, B2, and B3, respectively (Fig.Q^d) is a snapshot 
sequence of a Model Bl run). The general trends seen in Fig. 21 are similar to those seen in Fig. |21 The 
head-on cases tend to result in spherical remnants of decreasing mass with increasing v. The elongation of 
the remnants tends to increase with an increase in 6, up to a point. Of the three models note that Model Bl 
is the most similar to Model A. This is because Model Bl has the same amount of net angular momentum 
in the system since the spin components of the impactors cancel. Model B2 and Model B3, however, have 
smaller and larger net angular momentum in the system, respectively, than Model A or Bl. This is reflected 
in the number of runs with fast-rotating and/or elongated remnants f Table ITl)l . Models A and Bl have an 
intermediate number of runs with extreme P and/or e values compared with Model B2 or B3 (Model C is 
a special case discussed in the next section). 

Model Bl does differ from Model A in one important respect. As seen in Fig.Qfa), some of the remnants 
in this model (e.g., b — 0.30, v — 1.10; h = 0.60, v — 0.61) have unique asymmetries (i.e., broken eight-fold 
symmetry) . This is because before the encounter one of the bodies is spinning prograde while the other body 
is spinning retrograde with respect to the orbit. The prograde rotator has larger angular momentum with 
respect to the center of mass of the system than its retrograde counterpart, consequently, it suffers more 
mass loss. This is analogous to the resistance of retrograde rotators to tidal disruption (Richardson et al. 
1998). 

To summarize other quantitative results, 24% of the Model Bl runs resulted in net erosion, while this 
value was 29% for B2, and 40% for B3. The mixing statistics are generally similar to those for Model A, 
namely that larger disruption resulted in more mixing. As for ejecta statistics, again no more than about 
2% by mass remains in orbit around the remnant in all cases, while a somewhat larger percentage is destined 
to reaccrete (no more than 6%, except for a few cases where components of a future contact binary were 
on slow-return trajectories). The distribution of fragments (ni, n2, and n) followed similar trends to those 
of Model A. 



3.1.3 Model C: Unequal size, no spin 

In Model C we used two different-sized impactors witli no initial spin: one large sphere of 1357 particles 
and 1 km radius, and one small sphere of 717 particles and 0.46 km radius, keeping the bulk densities the 
same (2 g cm^"^) and the total number of particles similar to the previous models. Hence the larger sphere 
is ten times the mass of the smaller sphere and the particles in the two impactors are different sizes (the 
smaller body has smaller particles, to ensure adequate resolution). We caution that the difference in particle 
sizes implies different packing efficiencies (porosities) which may affect the outcome (c/. Section rO|l . Both 
impactors were equilibrated using the same process as before. Note that the parameter space investigated is 
different from the previous models, primarily for better sampling of the tidal regime (large b, small v). For 
this model, Wcrit = 2.9 m s^^. As for the previous models, the minimum v is twice the value corresponding 

to Voo = 0. 

In Fig. Efd) it is evident that most collisions result in net growth of the larger body (only 9 cases, or 
16%, result in net erosion i.e., remnants with less than 90% of the total mass of the system). None of the 
encounters resulted in critical dispersal and only the highest-speed cases resulted in violent disruption of the 
combined system (i.e., b ^ 0, v — 2.5). Also the remnants are all roughly spherical (the largest e = 0.26). 
Fig. ^e) shows a typical encounter: the small body is pulverized and occasionally planes off a chunk of the 
larger body (so ni is typically a few hundred in all runs except the most grazing, while n2 and n remain 
small, J; 10). Most of the smaller fragments escape the system, a tiny fraction (< 1% by mass) go into orbit 
around the remnant, while the remaining fragments return and blanket the remnant in the equatorial plane. 
The largest concentration of smaller particles is at the impact site. The rotation periods of the remnants in 
this model are typically long compared to those of the previous models due to the larger rotational inertia 
of the bigger impactor. Finally, there was little tidal interaction seen in any of the cases, suggesting even 
the minimum v was too large. Unfortunately, smaller v would result in stronger path deflections, making 
interpretation difficult. 

3.2 CoefRcient of Restitution Test 

The energy change in the center-of-mass frame of a system of two smooth, colliding spheres is given by (e.g., 
Araki and Tremaine 1986): 

AE = -^ti{l~el)vl (8) 

where Vn is the component of relative velocity normal to the mutual surfaces at the point of contact, and 
/I and e„ have the usual definitions. Hence as e„ 0, all the impact energy — less a geometric factor that 
depends on b — is dissipated. Although a collision between two rubble piles consists of many individual 
particle collisions, the dependence of AE on e„ suggests that will depend on e„ in a similar way, namely 
that smaller e„ implies larger Q^. 

A simple test bears this out. Table IHII and Fig. [3 summarize the effect of varying e„ for one of the Model 
A runs {b = 0.15, v = 2.00; cf. Fig. ^b)). The general trend is clear: as e„ decreases, the size of the largest 
remnant increases (note the large Afacc value for the e„ — 0.2 case, indicating that a big fragment is about 
to merge with the remnant, giving it the largest mass of all the runs). Runs with smaller e„ form discrete 



rubble piles out of the collision debris faster and more efficiently than those with larger e„. For e„ = 1, no 
rubble piles actually form. The strong dependence on e„ suggests that further study is needed to determine 
the value most representative of true rubble pile collisions. 

3.3 High-Resolution Models 

In order to test the degree to which particle resolution affects the collision outcome, we performed two 
high-resolution runs using parameters drawn from Model A. Each impactor for this test consisted of 4,995 
identical particles, more than 5 times the number used in the parameter space runs f Section 18.1(1 . The 
progenitor needed 1 CPU day to equilibrate using 2 processors on the SGI Origin 200. At equilibrium, 
the code was performing ~ 4 x 10^ collisions per step, with each step requiring ~ 140 s wallclock time. 
The enhanced packing efficiency of the high-resolution impactors plus their randomized orientations makes 
detailed comparison with the low-resolution runs difficult. However, we would expect the general trends to 
be similar (i.e., outcome class, etc.). Figure|HJa) shows snapshots shortly after the initial impact comparing 
the low- and high-resolution Model A runs with b = 0.30, v = 1.25. Figure El^b) shows post-reaccretion 
snapshots for b — 0.60, v — 0.61. These runs were chosen because they are moderately well separated in b-v 
space while still being representative of the complex intermediate-energy regime (c/. Fig. [21). The expense of 
these calculations precluded a more thorough sampling. 

Both high-resolution runs in this test show evolution similar to that of their low-resolution counterparts. 
In Fig.ini^a), the impactors mutually penetrate and lose most of their relative orbital energy (the bodies will 
eventually accrete into a single massive remnant). Note the presence of the "mass bridge" between the two 
bodies in both cases. The rotational phase and penetration distance differ somewhat, perhaps indicating that 
higher resolution (and hence lower porosity) gives rise to more efficient dissipation, perhaps by increasing the 
degrees of freedom. In Fig.|S{b), the final shape of the reaccreted body at low and high resolution is similar, 
but there is more structural detail in the high-resolution remnant, e.g., the depression on the upper surface. 
The remnant mass and rotation period are comparable at this instant: 0.985 and 4.2 h respectively at low 
resolution; 0.987 and 4.1 h at high resolution. We conclude that higher resolution may give insight into the 
more detailed aspects of reaccumulation, but low resolution is sufficient for a broad sampling of parameter 
space. 

4 DISCUSSION 

4.1 Critical Dispersal Threshold 

Despite the large number of runs carried out for this investigation, the data are still too sparse in each 
model to reliably derive a generalized expression for the retained mass (remnant plus accreting and orbiting 
material) as a function of b and v, i.e., 1 — M^sc = f{i>,v). However, we can solve for the critical contour 
/(&, v) — 0.5, which is well sampled by our choice of parameter space (for Model C, we solve for f{b, v) = 0.9, 
the point of net erosion for the larger impactor). Our method is to perform bi- linear interpolation of our 
b vs. V results onto a regular grid, root solve using Newton's method for the v value that gives a remnant 



mass of 0.5 at each grid line in b (we chose 20 hnes for smooth samphng), and fit the resulting values to a 
functional form. After some experimentation, we found the contour is best represented by a Gaussian: 



v\f=o.5 = Vi, = aexp 
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(9) 



where a, (3, 7, and S are parameters to be determined by non-linear least-squares fitting. 

Table Hvl gives the best-fit values of the Gaussian parameters along with their 1-a uncertainties for each 
of the parameter space models. Note that the fits are marginally consistent with (3 — 0, i.e., no b offset, 
except for Model C. The differences between the fits (except Model C) are slight, but they follow the trend 
mentioned in Section [3. 1.21 namely that Model B2 has a higher disruption threshold than Model B3, with 
Models A and Bl having intermediate thresholds. Model C has a broader distribution that is somewhat 
offset in &, but inspection of the other models shows similar trends for the 0.9 contour. 



4.2 Debris Size Distributions 

Combining all 275 runs of Models A, B, and C, we find the largest primary (remnant) mass is 0.999 (so 
there were no perfect mergers), the largest secondary mass is 0.498 (there was always some grazing mass 
exchange, at least in Models A and B), and the largest tertiary mass is 0.073. The smallest primary mass 
is 0.038. In Models A and B, the most common outcome was an even split in mass between the primary 
and secondary, since most runs at moderate to large b resulted in little to no mass exchange between the 
impactors. For b < 0.30 (0 < 17°), the normalized primary mass function is well approximated by a curve 
of the form n('m) cx 1/(1 — m?), m < 1 (Fig. 0I- 



4.3 Debris Spatial Distributions 

In cases where debris escapes the central remnant, the ejected material is invariably concentrated in a plane 
normal to the orbital (z = 0) plane, although in some cases material can be spread out in the orbital plane 
as two returning fragments coalesce. For our head-on collisions (6 = 0), the dispersal plane is normal to 
the X axis. For 6 > 0, the plane is initially normal to the impact angle 0, but rotational inertia from the 
orbital motion causes the dispersal plane to overshoot this value. As usual, initial spin may help or hinder 
this process (note for Model B3, (f> ^ — sin^^ b). 

Figure |S1 illustrates the anisotropic distribution of ejecta as projected to the orbital plane for 4 of the 5 
runs shown in Fig. ^ In the equal-mass cases the angular distribution is bi-modal, with the peaks roughly 
180° apart. For the plot representing Fig. ^b), the peaks are of unequal amplitude since the remnant is 
relatively small and displaced from the system center of mass. The one unequal-mass case is uni-modal, 
indicating that debris was scattered preferentially in one direction (roughly 30° measured counterclockwise 
from the x axis) , as seen in Fig. ^e) . 

Generally the z distributions are sharply peaked near the largest remnant but some particles end up 
many hundreds of km away. Recall that the tidal field of the Sun is not included in our simulations. If it 
were, these particles would be well outside the remnant's Hill sphere at 1 AU (for example): 
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Regardless, most of these particles would escape the remnant, even without the solar tides. 



4.4 Outcome Probability 

For a given impact parameter and speed distribution, the probability of a net accretional (as opposed to net 
erosional) outcome can be estimated from Eq. J^. Suppose we set b — 0.7, which corresponds to ip — 45°, 
the most probable impact angle for randomly flying projectiles striking a spherical target (Love and Ahrens 
1996). A monodispersive population of bodies with a Maxwellian distribution of speeds whose rms equals 
the escape speed Ve from a particle's surface has the following normalized distribution function in relative 
speed (e.g., Binney and Tremaine 1987, Problem 7-3): 

9i'")dv = 1_ exp (-^^ v^dv. (11) 



2^vi ^ V 41-2 ; 

The probability of a net accretional impact for hyperbolic encounters with 6 = 0.7 is then: 

q(v)dv 

P[f{b = 0.7, v) > 0.5] = 7^'^ ' , (12) 

where Vo is the initial speed corresponding to I'oo = from Eq. Q and v-k is obtained from Eq. (j^J. For 
Model A, we have — 0.51, Vo — 0.22, and = 0.81. Solving Eq. H12I) numerically we find the probability 
of an accretional impact in this case is 26%. The probability of erosion is 1 — P{f > 0.5) = 74%. 

The full accretional cross section is obtained by integrating Eq. (|12|l over all impact parameters. The net 
accretion probability is then the ratio of this value to the geometrical cross section: 

P[f{b,v) > 0.5] = — T— ^ / bdh '"^'^ , (13) 

nj^bdb Jo l^^,^g{v)dv 

where the dependence of Vo and v^, on b has been made explicit. Solving this equation we find the accretion 
probability increases to 36% only, since head-on collisions are relatively rare. Such a low value implies that 
this population of rubble piles would not go on to form planets but would instead grind itself down to dust. 
If rubble piles were common during the early stages of planet formation then perhaps collisions were more 
dissipative than modeled here. Alternatively, the accretion probability may be enhanced when there is a 
distribution of masses, a possibility that can only be tested with more simulations using impactors of varying 
size. 

4.5 Comparison with Previous Work (Gravity Regime) 

Using the fits to Eq. ^ we can estimate the value of (recall this is for e„ — 0.8; further runs are needed to 
determine the dependence on dissipation). Restricting ourselves to Model A with 6 = 0, we find ~ 1.9 J 
kg~^. This lies very close to the Holsapple (1994) and gravitational binding energy curves described in Love 
and Ahrens (1996; see in particular their Eq. (2) and Fig. 7). It lies well off the extrapolation of their SPH 
results. In their paper they suggest that the discrepancy between their results and analytic or experimental 
results may arise from: 1) the local rather than global deposition of impact energy at the surface of the 
target; 2) the difference between the role of gravity in self-compression and ejecta retention; 3) the finite size 
of the projectile. The present work however is similar to Love and Ahrens (1996) in all these respects, which 



suggests the difference may be attributable solely to the adopted equation of state (an incompressible fluid 
in our case, compared with the Tillotson equation of state for granite in theirs). Note that the SPH curve 
plotted in Fig. 7 of Love and Ahrens (1996) was for an impact angle oi <f) = 45°, but this would amount to 
less than an order of magnitude difference in Q^. 

If the outcome truly depends solely on the gravitational binding energy (ignoring the effect of dissipation 
for now as this requires further study), then we would expect Q*j-, oc M/R oc R'^ oc M^/^. Prom our Model 
A point we can estimate the constant of proportionality: Q*j-, 1.2 x 10~^i?^ 2.9 x 10~^M^/^. Further 
models with different M are needed to confirm this result (our Model C case failed to sample the critical 
dispersal regime, so we cannot use it here). 

Watanabe and Miyama (1992) found Mesc oc for their low-speed, head-on SPH models (see Eq. (3.5.1) 
in their paper). We find a similar trend. For the 6 = outcomes of Model A, a least-squares fit to the form 

Mesc = avf^ (14) 

yields a = 0.06 ± 0.02 and /3 = 3.2 ± 0.1. Evidently this relation must break down at large v, otherwise Mesc 
would exceed unity. Indeed our only significant outlier is for our highest v value (2.50), with Mesc in this 
case 20% below the curve. 

5 CONCLUSIONS 

In summary, we have conducted a series of numerical simulations to create a partial map of the parameter 
space of rubble pile collisions at low impact speeds. The general trends can be summarized as follows: 
1) larger impact angles result in more elongated, faster-spinning remnants; 2) larger impact speeds result 
in greater mass loss and increased mixing of the remnant; 3) initial impactor spin can increase or reduce 
the rotation period and elongation of the remnant. It is also possible to create asymmetric shapes if the 
impactors have oppositely oriented spins. These general trends are directly related to the total energy and 
angular momentum of the system. In cases where one impactor is significantly larger than the other (Model 
C), the smaller body generally disrupts completely on impact, sometimes removing a modest fraction of the 
surface of the target body and sometimes redepositing material along the remnant's z = equator. 

We have been able to generate a wide variety of remnant shapes, including spheroids, ellipsoids, contact 
binaries (peanut shaped and S shaped), and shapes with broken eight-fold symmetry. It proved considerably 
difficult to get a significant amount of material to orbit the remnant; most debris (98%) either accreted 
onto the remnant or escaped from the system. We found no detached binaries of significant size, but ~ 
10% of the remnants in Model A and B arc contact binaries. The coefficient of restitution appears to play 
a more important role in collisions than in tidal disruption and can strongly affect the number and size of 
post-impact rubble-pile fragments. Increased particle resolution (or reduced porosity) appears to augment 
dissipation and give rise to more complex shapes, but the effects are modest over a factor of 5 in particle 
number. 

We found that the impact speed needed for critical dispersal is well represented by a Gaussian function of 
impact parameter. Given a velocity distribution it is possible to estimate the probability of cither impactor 
gaining or losing mass as a result of the collision. At low impact angles with equal-size impactors the remnant 



mass function is roughly proportional to 1/(1 — m^). Secondary and tertiary masses are typically finite but 
small, except for near-grazing encounters. Most material is ejected anisotropically in a plane perpendicular 
to the axis of the initial motion and in some cases the debris can coalesce into smaller rubble piles. 

We found Q}} 2 J kg~^ for the head-on collisions in Model A and that Mggc oc v^'"^. The former result 
is in rough agreement with the theoretical gravity-regime model of Holsapple (1994). The latter relation 
agrees with Watanabe and Miyama (1992). We find that km-sized rubble piles in general are much easier 
to disperse than extrapolation of the strengthless granite models of Love and Ahrens (1996) would suggest. 
This may be due in part to our conservative choice for e„ (0.8). Although more work needs to be done, we 
believe our simulations may provide a numerical basis for parameterizing collisions during the early stage of 
planet formation, when the planetesimals are dynamically cool and the dominant sizes are still $ 10 km. 

5.1 Future Work 

In this study wc were restricted to investigating the dependence of collision outcome primarily on impact 
parameter and impact speed. We also examined a few spin combinations, a model with unequal masses, 
and a single run with various values of the restitution coefficient. But the parameter space is truly vast. 
Naturally we would like to test more values for the parameters we have already investigated, particularly the 
coefficient of restitution and the dependence on porosity. We also need a finer grid at small speed and near- 
grazing separation to fully investigate the tidal regime (this would also provide better data for comparison 
with stellar system collision models, e.g., Davies et al. 1991). However there are many other new parameters 
to explore. We would like to test the effect of changing the spin-axis orientations (beyond pure prograde 
or retrograde). We suspect this would result in even more unusual shapes. Non-spherical impactors with 
a variety of sizes would improve realism and provide better estimates of Q^. A spectrum of particle sizes 
could alter the effective dissipation as smaller particles fill the voids between larger ones. Adding surface 
friction could lead to steeper slopes and the possibility of simulating crater formation in large targets. We 
plan to add a simple model for compaction to allow higher impact speeds and compare with the porous 
models of Housen and Holsapple (1999). We would like to track the movement of particles near the cores of 
our nibble piles and compare with the surface particles to study "scrambling" in a single rubble pile. There 
are so many possibilities that likely the only practical approach would be to randomly sample points in this 
vast parameter space to get a feel for the overall trends and then concentrate on the most interesting aspects 
in detail. We will carry out such work in the future. 
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Table I: Summary of Model A results (Section 13.1. l|l 
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Table II: Comparison of runs with extreme P and e values (Section 
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No. with P < 5 h 


No. with e > 0.35 
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11 
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10 
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10 
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15 


10 
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Table III: Effect of varying dissipation in Model A run b — 0.15, v — 2.00 (Section l3.2|l 
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Table IV: Parameter fits to Eq. (jg)) for Mrem = 0.5 contours, b < 1 
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Figure Captions 

Figure U Snapshots of rubble pile collisions from representative runs as seen in the center-of-mass frame. 
The models and runs are: (a) Model A, 6 = 0.00, v = 1.00; (b) Model A, b ^ 0.15, v ^ 2.00; (c) Model 
A, 6 = 0.90, V = 0.52; (d) Model Bl, b = 0.30, v = 1.10; and (e) Model C, 5 = 0.50, v = 1.25. The 
arrow of time is to the right. The interval between frames is not regular: the snapshots were chosen 
to highlight distinct stages in the evolution of each run. In run (b), the final two frames have been 
brightened for clarity. 

Figure [2| Projected shape of the largest remnant at the end of each Model A run as a function of b and 
V. At this scale each grid square measures 4 km on a side. Solid inner squares indicate remnants that 
retain at least 90% of the system mass; dashed squares indicate remnants with at least 50%. Critical 
dispersal generally corresponds to the transition from solid to dashed, although in some cases a sizeable 
fragment may be about to accrete with the remnant. Table Ogives additional data for this model. 

Figure (S) Illustration of the spin sense for the Model B impactors. In Model Bl, the impactors have 
opposite spin; in B2 they have the same spin, oppositely aligned with the orbital angular momentum; 
in B3 the spins are aligned with the orbital angular momentum. 

Figure 131 Remnant shapes for the remaining parameter space models. The model is indicated in the top 
left of each plot. Compare with Fig. [5] 

Figure |5l Snapshots showing the effect of varying e„ for the Model A run with b — 0.15, v — 2.00 (c/. 
Fig. n^b)). Each snapshot was taken about 6.5 h after impact. The chosen e„ value and camera 
zoom-out factor are shown in the top left of each frame. For clarity, no color or shading distinction 
is made between the particles of the original impactors, and the e„ — 0.9 and e„ = 1.0 frames have 
been brightened. From these snapshots and the statistics in Table Hill it can be seen that rubble pile 
formation favors smaller e„ values. The differences at the extremes are dramatic. 

Figure |6j Comparison of two Model A runs performed at low resolution (955 particles per rubble pile; 
left column) and high resolution (4,995 particles per rubble pile; right column). The nm parameters 
are: (a) b = 0.30, v = 1.25; (b) b — 0.60, v — 0.61. The evolution is similar in both cases, with 
differences attributable to packing efficiency, initial orientation, and possibly enhanced dissipation at 
higher resolution. 

Figure [7| Primary (solid line) , secondary (dotted) , and tertiary (short dashed) mass fractions of every 
Model A and B nm with b < 0.30, sorted by primary mass. The long-dashed line was obtained by 
integrating a rough fit to the primary mass function of the form 1/(1 — m^). 

Figure |8j Debris dispersal patterns in the initial orbital plane relative to the largest remnant for the runs 
labeled (a), (b), (d), and (e) in Fig. ^ The 9 histograms are binned in 5° increments. Only particles 
with projected distances in the z plane exceeding twice the remnant radius were included. 
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